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Abstract 

The instantaneous state of a neural network consists of both the degree 
of excitation of each neuron the network is composed of and positions of 
impulses in communication lines between the neurons. In neurophysiolog- 
ical experiments, the neuronal firing moments are registered, but not the 
state of communication lines. But future spiking moments depend essen- 
tially on the past positions of impulses in the lines. This suggests, that 
the sequence of intervals between firing moments (inter-spike intervals, 
ISIs) in the network could be non-Markovian. 

In this paper, we address this question for a simplest possible neu- 
ral "net", namely, a single inhibitory neuron with delayed feedback. The 
neuron receives excitatory input from the driving Poisson stream and in- 
hibitory impulses from its own output through the feedback line. We 
obtain analytic expressions for conditional probability density P{t„+i \ 
t„, . . . ,ti,to), which gives the probability to get an output ISI of du- 
ration tn+i provided the previous (n + 1) output ISIs had durations 
tn, ■ ■ ■ ,ti,to. It is proven exactly, that P{tn+i \ tn, . . . , ii, io) does not 
reduce to P(t„+i | t„, . . . for any n > 0. This means that the output 
ISIs stream cannot be represented as a Markov chain of any finite order. 

1 Introduction 

In a neural network, the main component parts are neurons and inter-neuronal 
communication lines - axons \Xi- These same units are the main ones in most 
types of artificial neural networks [S]- If so, then the instantaneous dynamical 
state of a network must include dynamical states of all the neurons and com- 
munication lines the network is composed of. The state of a neuron can be 
described as its degree of excitation. The state of a line consists of information 
of whether the line is empty or conducts an impulse. If it does conduct, then 
the state of the line can be described by the amount of time which is required 
for the impulse to reach the end of the line (time to live). 

In neurophysiological experiments, the triggering (spiking, firing) moments 
of individual neurons but not the states of communication lines are registered. 
The sequence of intervals between the consecutive moments (inter-spike inter- 
vals, ISIs) is frequently considered as a renewal [3] or Markovian [4| stochastic 
process. For a renewal process, the consecutive ISIs are mutually statistically 
independent. Moreover, all statistical characteristics of a spike train must be 
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derivable from the single-ISI probability distribution. Additionally, those char- 
acteristics must be the same for a shuffled spike train, obtained by randomly 
reordering the ISIs, since shuffling does not change the single-ISI probability 
distribution. On the other hand, the experimentally obtained spike trains in 
auditory ^ and visual [6] sensory systems does not support the ISIs' mutual in- 
dependence. This is revealed by calculating the correlation coefficient between 
the adjacent ISIs, which appeared to be nonzero for the experimental spike 
trains, while it must be zero for any renewal process. Also, such characteristics 
as Fano factor curve and firing rate distribution calculated for shuffled spike 
trains differ qualitatively from those obtained for the intact ones. These obser- 
vations can be associated with memory effects in the ISI sequence which arise 
from an underlying non-renewal process. Recently [7], such a possibility was an- 
alyzed for weakly electric fish electrosensory afferents using high-order interval 
analysis, count analysis, and Markov-order analysis. The authors conclude that 
the experimental evidence cannot reject the null hypothesis that the underlying 
Markov chain model is of order m or higher, or maybe non-Markovian. The 
limited data sets used in allow to establish a lower bound for m as m > 7 
for some fibers. 

What could be possible sources of such non-renewal, or even non-Markovian, 
behavior in real neural network? First, this behavior could be inherited from 
non-renewal (non-Markovian) character of the input signal. Second, intrinsic 
neuronal properties, such as adaptation, could be responsible. Finally, as we 
show here, the presence of delayed feedback interconnections itself could be the 
possible source of the non-Markovian behavior. 

The non-Markovian behavior of the ISI sequence from neuron in a network 
with delayed interconnections is not surprising. Indeed, the information about 
which neurons are spiking/silent at any given moment of time leaves unknown 
the position of impulses in the interconnection lines at that moment. And it is 
the previous firing moments which determine the states of interconnection lines, 
which in turn determine the next firing moments. Therefore, information about 
the previous neuronal firing moments could improve our predicting ability as 
regards the next firing moments. 

In this paper, we consider a simplest neural "net" , namely, a single neuron 
with delayed feedback, which is driven with Poisson process. As neuronal model 
we take binding neuron as it allows rigorous mathematical treatment. We study 
the ISI output stream of this system and prove that it cannot be presented as 
Markovian chain of any finite order. This suggests that activity of any network 
with delayed interconnections, if presented in terms of neuronal firing moments, 
should be non-Markovian as well. 

2 The object under consideration 
2.1 Binding neuron model 

The understanding of mechanisms of higher brain functions expects a contin- 
uous reduction from higher activities to lower ones, eventually, to activities in 
individual neurons, expressed in terms of membrane potentials and ionic cur- 
rents. But the description of the higher brain functions in terms of potentials 
and currents in parts of individual neurons would be difficult, similarly as it 
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Figure 1: Binding neuron with feedback hne under Poisson stimulation. Multi- 
ple input lines with Poisson streams are joined into a single one here. A is the 
delay duration in the feedback line. 



would be difficult to describe execution of computer programs by a CPU in 
terms of KirhgofF's laws. In this connection, it would be helpful to abstract 
from the rules by which a neuron changes its membrane potentials to rules by 
which the input impulse signals are processed in the neuron and determine its 
output firing activity. The coincidence detector, and temporal integrator are 
the examples of such an abstraction, see discussion in [S] . 

One more abstraction, the binding neuron (BN) model, is proposed as sig- 
nal processing unit |10) . which can operate either as coincidence detector, or 
temporal integrator, depending on quantitative characteristics of stimulation 
applied. This conforms with behavior of real neurons, see, e.g. ]T1\ [T^. The 
BN model describes functioning of a neuron in terms of discrete events, which 
are input and output impulses, and degree of temporal coherence between the 
input events, see [5] for detailed description. Mathematically, this model can 
be realized as follows. We expect that all input impulses in all input lines are 
identical. Each input impulse is stored in the BN for a fixed time, t. The r is 
similar to the tolerance interval discussed in |13j . All input lines are excitatory. 
The neuron fires an output impulse if the number of stored impulses, E, is equal 
or higher than the threshold value, Nq. After that, BN clears its memory and 
is ready to receive fresh inputs. That is, every input impulse either disappears 
contributing to a triggering event, or it is lost after spending r units of time in 
the neuron's internal memory. The latter represents leakage. Here, the leakage 
is abrupt, while in more traditional models it is gradual. 

The BN model is not general, but somewhat inspired by neurons as integra- 
tors up to a threshold. Its name is suggested by binding of features/events in 
large-scale neuronal circuits [141 115[ 116) . Its operational simplicity is provided 
by the fact that each input impulse traces entirely disappear after finite time 
r. This is in the contrast to more familiar models where the traces (excita- 
tory postsynaptic potentials, EPSP) decay exponentially. E. g., in the leaky 
integrate-and-fire model, EPSP is mimicked as pure exponential function the 
traces of which can disappear completely only after triggering. In the BN model, 
the EPSP is mimicked as box function of width/duration r and the traces are 
stored in the neuron no longer than r units of time. 

Further, we expect that input stream in each input line is the Poisson one 
with some intensity A^. In this case, all input lines can be collapsed into a single 
one delivering Poisson stream of intensity A = Ai, see Figure [TJ 

For analytic derivation, we use BN with iVo = 2 in order to keep mathemat- 
ical expressions shorter. It seems, that cases with higher thresholds might be 
considered with the same approach, but even A'o = 3 without feedback requires 
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additional combinatorial efforts, see [T^. Therefore, cases of higher threshold 
are tested here only numerically. 

As regards real biological neurons, the number of synaptic impulses in the 
internal memory which is necessary to trigger a neuron, varies from one [18) . 
through fifty [II], to 60-180 f20l, and 100-300 [21]. 

2.2 Feedback line action 

In real neuronal systems, a neuron can form synapses from its axonal branch 
to its own dendritic tree [HHKHHnilinillZlllHlIlSI- Synapses of this type 
are called autapses. Some of the neurons forming autapses are known to be 
inhibitory, see [2H [27J [55] for experimental evidence. As a result, the neuron 
stimulates itself obtaining an inhibitory impulse through an autapse after each 
firing with some propagation delay. We model this situation assuming that 
output impulses of BN are fed back into BN's input with delay A. This gives 
the inhibitory BN with delayed feedback model. Figure [TJ 

The inhibitory action of feedback impulses is modeled in the following way. 
When the inhibitory impulse reaches BN, it annihilates all excitatory impulses 
already present in the BN's memory, similarly as the Cl-type inhibition shunts 
depolarization of excitable membrane, see |31j . If at the moment of inhibitory 
impulse arrival, the neuron is empty, then the impulse disappears without any 
action, similarly as Cl-type inhibition does not affect membrane's voltage in its 
resting state. Such inhibition is "fast" in that sense, that the inhibitory impulses 
act instantaneously and are not remembered by neuron. This simple behavior 
is approved by relatively fast kinetics of the chloride inhibitory postsynaptic 
currents [35]. 

The feedback line either keeps one impulse, or keeps no impulses and cannot 
convey two or more impulses at the same time. Biological correlates supporting 
to an extent this assumption could be a prolonged refractory time and/or short- 
term synaptic depression. The latter can have the recovery time up to 20 s 
[30| . If the feedback line is empty at the moment of firing, the output impulse 
enters the line, and after time interval equal A reaches the BN's input. If the 
line already keeps one impulse at the moment of firing, the just fired impulse 
ignores the line. 

This means, that at the beginning of an output ISI the feedback line is 
never empty. In order to describe the state of the feedback line, we introduce 
the stochastic variable s, s g]0; A], which gives the time to live of the impulse 
in the feedback line, see Fig [T] Hereinafter, we will use the values of s just at 
the moments of output ISI beginnings (just after firings). 

We assume, that time delay A of impulse in the feedback line is smaller than 
the BN's memory duration, r: 

A < r. (1) 

It allows to make analytic expressions shorter. Also, the assumption ([T]) is 
consistent with the case of direct feedback, not mediated by other neurons. 

3 Statement of the problem 

The input stream of impulses, which drives neuronal activity is the Poisson 
stream. It is stochastic, therefore, the output activity of our system requires 
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probabilistic description in spite of the fact that both the BN and the feed- 
back line action mechanisms are deterministic. We treat the oiitput stream of 
inhibitory BN with delayed feedback as the stationary procesqj- In order to 
describe its statistics, we introduce the following basic functions: 

• the joint probability density P{t„i,tm~i, ■ ■ ■ ,to) for (m + 1) successive 
output ISI durations, to is the first one. 

• the conditional probability density P{tm \ tm-i, • • ■ , ^o) for output ISI du- 
rations; P{tm I tm-l^ ■ ■ ■ , to)dtm givcs the probability to obtain an output 
ISI of duration between t„i and tm + dim provided the previous m ISIs 
had durations t„i-i,t„i-2, ■ ■ ■ ,to, respectively. 

Definition 1. The sequence of random variables {tj}, taking values in Vl, is 
called the Markov chain of the order n > 0, if 

Vm>nVtogn . . . Vf^go P{tm \ tm-l, • ■ • , ^o) = P{tm \ tm-l, • ■ • , t„i-n), 

and this equation does not hold for any n' < n (see e.g. \8S^). In the case of 
ISIs one reads fl = M"*" . 

In particular, taking m — n + 1, we have the necessary condition 

P{tn+1 I tji, ■ ■ . ,ti,to) = P{tn+l I tn, . . . ,ti), 

gR+, i = 0,...,n + l, (2) 

required for the stochastic process {tj} of ISIs to be the n-order Markov chain. 
Our purpose in this paper is to prove the following theorem. 

Theorem 1. The output ISIs stream of inhibitory BN with delayed feedback 
under Poisson stimulation cannot be represented as a Markov chain of any finite 
order. 

4 Main calculations 

This section with Appendices contains the required proof of Theorem [1] 
4.1 Proof outline 

In order to prove the Theorem [U we are going to show analytically, that the 
equality ^ does not hold for any finite value of n. Namely, we will derive 
the exact analytic expression for the conditional probability density P(t„+i | 
tn, . ■ ■ ,ti, to) and show, that it depends on to for any finite number n. 
For this purpose we introduce the stream ts of events (t, s) 

ts — {. . . , (ti, 5^), . . . }, 

where Si is the time to live of the impulse in the feedback line at the moment, 

when the ISI ti starts. We consider the joint probability density P{tn+i, s„+i; i„, s„; . . . ; to, so) 

^ The stationarity of the output stream results both from the stationarity of the input 
one and from the absence of time-dependent parameters in the BN model, see Section 12.11 
In order to ensure stationarity, we also expect that system is considered after initial period 
sufficient to forget the initial conditions. 
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for realization of (n + 2) successive events (i, s), and the corresponding condi- 
tional probability density P{tn+ij Sn+i \ tn, s„; . . . ; toi so) for these events. 

Then, we proof the following lemma, which will be used in our calculations. 

Lemma 1. Stream ts is the 1-st order Markovian: 



V„>oVto>oVs„g]0;A] ■ • ■ Vt„_|_i>oVs„_^^g]0;A] 

P{tn+1, Sn+i I tn, S„; . . . ; tg, Sq) — P{tn+1, Sn+1 | tn, Sn), (3) 

where {to, . . . ,tn+i} is the set of successive ISIs, and {sq, . . . , Sn+i} are the 
corresponding times to live. 

See Appendix El for the proof. 

Then, in order to find the conditional probability density P{tn+i \ tn, ■ ■ ■ ,ti,to), 
we perform the following steps: 

• Step 1. Use the property Q for calculating joint probability density of 
events {t, s): 

P{tn+1: Sn+l 

P{tn+l,Sn+l I tn, S„) ■ ■ • ^'(^1,51 | ^0, SQ)P{tQ, Sq), (4) 

where P{t,s) and P(tn,Sn \ in-i,s„_i) denote the stationary probabil- 
ity density and conditional probability density (transition probability) for 
events {t,s). 

• Step 2. Represent P{tn+i,tn, . . . ,to) as marginal probability by integra- 
tion over variables s,;, j = 0, 1, . . . , n -I- 1: 

P{tn+1, tn, ■ ■ ■ , to) = 

/ dso / dsi . . . / dSn+lP{tn+l,Sn+i;tn,Sn;.-.;to,So). (5) 

Jo Jo Jo 

• Step 3. Use the definition of conditional probability density: 

P{tn+1 \tnT ■ ■ Tti,to) ^ — ''p^^' ■ (^) 

Taking into account the Steps 1 and 2, one derives for the joint probability 
density 

P{tn+1, tn, ■ ■ ■ , to) — 

pA fA "+1 

/ dsQ.../ dsn+iP{to,so) P(tfe,Sfc I tfc_i,Sfe_i). (7) 
Jo Jo 

In the next sections, we are going to find the exact analytic expressions 
for probability densities P{t,s) and P{tk,Sk \ tk-i, Sk-i), and perform the in- 
tegration in (O. Then we will apply the Step 3, above, to find expressions 
for the conditional probability densities P(tn+i \ tn, ■ ■ ■ ,to). It appears, that 
P{tn+i I tn, . . . ,to) is a function with jump discontinuities. In order to prove 
that the equality ^ does not hold for any n > 0, we analyze the positions of 
that jump discontinuities only. 
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4.2 Probability density P{t, s) for events (t, s) 

The probability density P{t, s) can be derived as the product 

P{t,s)^Fit\s)fis), (8) 

where /(s) denotes the stationary probabihty density for time to Hve of the 
impulse in the feedback line at the moment of an output ISI beginning, F{t \ s) 
denotes conditional probability density for ISI duration provided the time to live 
of the impulse in the feedback line equals s at the moment of this ISI beginning. 
Exact expressions for both /(s) and F{t \ s) for inhibitory neuron with delayed 
feedback are calculated in [33]. Particularly, 

F{t I s) = i (9) 
[(1 + As) e-^"pO(<-s), t>s, 

were P^{t), t > 0, denotes an output ISI probability density for BN without 
feedback, which was obtained in 17, Eq. (3)]. Explicit expressions for P'^{t) 
are different for different domains of t. For example, 

P°{t) = X'^t e^^\ t€]0;T]. (10) 

It is proven in [T^, that P'^{t) is a continuous function for whole range of ISI 
durations: t 6]0;oo[. 

It is essential for further study, that F{t \ s) considered as function of t has 
a jump discontinuity ai t = s. Indeed, using ^ and dT(I| . one obtains 

lim F{t I s) = A^s e"^" > 0, s e]0; A], 

t—>-s — 

lim F(t I s) = 0. 

The presence of jump in F{t \ s) at t — s can be explained as follows. Ac- 
cording to the definition of F{t \ s), the inhibitory impulse from the feedback 
line arrives s seconds later than the ISI t starts. After the inhibitory impulse 
arrival, it is guaranteed, that the BN is empty. To trigger the BN just after 
that moment, it is necessary to get two impulses from the input stream within 
infinitesimally small time interval. This event has infinitesimally small proba- 
bility for the Poisson process (as well as for any other point process). That is 
why, the value of probability density F{t \ s) drops to zero at i = s + and 
F(t I s) experiences discontinuity at t — s. 

It is important to emphasize, that F{t | s) is a continuous function elsewhere 
except of the point t — s, where it has strictly positive jump. This follows from 
(jni) and from the continuity of P'^{t). The continuity of F{t \ s) at t g]0; s[ and 
t e]s; oo[, and its jump at t = s will be used later. 

We also need an expression for /(s), which is 

4e2AA 

/(s)=a-<5(s-A) + <?(s), where a= (3 ^ 2AA)e2AA + i ' (^1) 

where 5{-) - is the Dirac delta-function, g{s) - is a regular function, which 
vanishes out of interval s g]0; A] (see [Ml Eq. (15)] for the exact expression). 
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Figure 2: Left: output ISI probability density P{t) reproduced from [Ml, Fig. 2]; 
Right: probability density /(s) for times to live of the impulse in the feedback 
line. Here r = 10 ms, A 8 ms, A 150 s^^, A^o==2. 



the a gives the probability to obtain the impulse in the feedback line with time 
to live equal A at the beginning of an arbitrary output ISI, A — is the input 
Poisson stream intensity. 

Let us explain the presence of Dirac (5-function type singularity in f{s). The 
probability to have time to live, s, exactly equal A at the moment of an output 
ISI beginning is not infinitesimally small. Every time, when the line is free at 
the moment of an output ISI beginning, the impulse enters the line and has 
time to live equal A. For the line to be free from impulses at the moment of 
triggering, it is necessary that t > s for the previous ISI. The set of realizations 
of the input Poisson process, each realization satisfying t > s, has non-zero 
probability a, see (|lip. and this gives the (5-function at s = A in the probability 
density /(s). 

The output ISI probability density P{t) for inhibitory neuron with delayed 
feedback can be obtained as the result of integration of ([S]) : 

Pit)= F{t\s)fis)ds. (12) 
Jo 

Discontinuity of F(t\s) at t — s and (5-function type singularity at s = A in /(s) 
result in discontinuity of P{t) at t = A. 

Examples of P{t) and /(s) graphs can be found in Fig. m 



4.3 Conditional probability density P{tk,Sk \ tk-i,Sk-i) 

Here we find the conditional probability density P{tk, Sk \ ife-i, Sk-i) for events 
{tk, Sk), which determines the probability to obtain the event (tk, Sk), with pre- 
cision dtkdsk, provided the previous event was {tk-i, Sk-i)- By definition of 
conditional probabilities, the probability density wanted can be represented as 
the following product 

P{tk,Sk I ife_i,Sfc_i) = F{tk I Sk,tk-i,Sk-i)f{sk I t/c-i, Sfc-i), (13) 

where F{tk \ Sk, ife-i, Sk-i) denotes conditional probability density for ISI dura- 
tion, tk , provided i) this ISI started with lifetime of impulse in the feedback line 
equal to s^, and ii) previous (t, s)-event was (ife_i,Sfc_i); the f{sk \ tk-i, Sk^i) 
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denotes conditional probability density for times to live of impulse in the feed- 
back line under condition ii). It is obvious, that 

F{tk I Sk,tk-uSk^i)^F{tk I Sfe), (14) 

because with Sk being known, the previous event (tk-i, Sk-i) does not add any 
information, useful to predict tk (compare with the proof of Lemma[Tl Appendix 

In order to find the probability density /(sfc | tk-i, Sk-i), let us consider 
various possible relations between tk-i and Sfc^i. If ife-i > Sk-i, the line will 
have time to get free from the impulse during the ISI tk~i. That is why at the 
beginning of the ISI tk, an output spike will enter the line and will have time to 
live Sfc = A with probability 1. Therefore, the probability density contains the 
corresponding (5-function: 

/(sfc I ifc-i,Sfc_i) = 5(sfc - A), tk-i>Sk-i. (15) 

If tk-i < Sk-i, than the ISI tk-i ends before the impulse leaves the feedback 
line. Therefore, at the beginning of the t^, the line still keeps the same impulse 
as at the beginning of tk-i- This impulse has time to live being equal to Sk — 

Sfc_i - so 

f{sk\tk-i,Sk-i)^5{sk-Sk-i+tk-i), ifc-i < Sfc-i- (16) 

Taking all together, for the conditional probability density P{tk, s/c | tk-i, Sk-i) 
one obtains 

P{tk,Sk I tfc_i,Sfc-i) F{tk I Sfe)(5(sfc - A), tk-i > Sfe_i, 

= F{tk I Sk)S{sk ~ Sfc-i + tk-i), tk-i < Sfc-i, (17) 

where exact expression for F{t \ s) is given in 



4.4 Joint probability density P(t„+i, . . . ,io) 

In this section, we are going to find the exact analytic expression for the joint 
probability density P(i„+i, . . . ,to) at the following domain 

Di = i^{to,...,tn,t„+i) |^t,<A|. (18) 

Notice, that coordinate tn+i is not included to the condition here. The set of 
(n + 2) successive ISI durations to, ... , tn,tn+i has non-zero probability, pA > 0, 
to fall into the domain p^ . Indeed, BN with threshold Nq — 2 requires 2{n + l) 
input impulses within time window ]0; A[ to be triggered (n + 1) times within 
this window (condition ([1]) ensures that no one input impulse will be lost). BN 
receives excitatory impulses from the Poisson stream and inhibitory impulses 
from the feedback line. But no more than one impulse from the line may have 
time to reach BN's input during time interval less than A. Therefore, if as much 
as (2n 4- 3) input impulses are received from the Poisson stream during the time 
interval ]0; A[, the inequality (|18p holds for sure, no matter was an impulse from 
the feedback line involved, or not. Therefore, pA > p{2n -I- 3, A) > 0, where 
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p{i, A) gives the probability to obtain i impulses from the Poisson stream during 
time interval A [33]: p{i,A) = e'^'^{XAy 

For a fixed (n + 2)-tuple (<o, ■ • • , ^n, Ui+i) G Di, let us split the integration 
domain for sq in ([7|) in the following way: 

]0; A] =]0; io]U]to; io + hpjto + ii; io + <i + ^2] U ■ • • U]to + ii + • • • + t„; A], 

or 

/ dso = / dso + 2^ / dso + / dso, 

and introduce the following notations: 

/ dso / dsi . . . / ds„+iP(io,so) J]^ P(ifc, | tfe-i, 

i = 0,l,2,...,n, (19) 

dso / dsi... / dsn+iP{to,so) J]^ P(tfe, | Sfe_i), (20) 



A;=l 



3 = 



where we assume, that X^jLji — for Ji > i2- 

According to dTj), (fT9)) and (pO)) . the probability density P{tn+i, ■ ■ ■ ,to) can 
be obtained as 

n+l 

P(t„+i,...,io) = ^/.. (21) 

i=0 

Substituting P{tQ,so) and P{tk,Sk \ tk-i, Sk-i) from expressions ([8]) and (flT)) 
to ([TO]) and (pO)) and performing integration over variables si,...,s„+i, one 
obtains 

n+l k-1 ^^=°*^ i fc-1 

- n ^(^'^ I A - E / n ^(^'^ I ^0 - E ^.■)5(so)ciso, 

z = 0,l,2,...,n. (22) 

n+l n+l k-1 

^„ f k=0 j=Q k=0 j=Q 

where F{t \ s) and g{s) were defined in ^ and pTI) (see Appendix [B] for the 
details of integration). 



In+l — 
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Taking into account (PT|) . ([^ and (^5]) . one obtains the following expression 
for the joint probability density for output ISI durations: 



n+1 



^'(^n+1, • ■ • , ^o) — E 

i=0 

n n+1 ,x J. „ t ,^ - J. 

= E n ^(^^ I ^ - E ^^■) / 5(so) n I ■'^0 - E 

i=0 fe=i+l j=i+i fc=0 J=0 

/'^ n+1 fc-1 n+1 k-1 

5(so) n I ^0 - E ^j)'i*o + « n ^(^fe I ^ - E ^j)' 

7, n n 7 n ^' n 



i=0 

V' t ■ 

n n+1 fc-1 ^ i fc-1 



1=0 



fe=o i=o fe=o j=o 



Y,t^<^^ n = 0,l,..., (24) 

i=0 

where we assume, that X^jLji — ^^^^ IIjLji — 1 ^o^' Ji > J2- 

The expression gives the joint probability density P(i„+i, . . . , to) for 
consecutive ISI durations at the domain Di for an arbitrary n. Therefore, 
the conditional probability density P{tn+i \ tn, . . . ,to) at Di can be obtained 
readily, see equation 

4.5 Discontinuities in P{tn+i, . . . ,to) 

In this section, we will answer two following questions; i) does the P(t„+i, . . . ,to) 
contain discontinuities at -Di? and ii) if it does, what are the positions of that 
discontinuities? 

In order to ascertain the continuity of expression, defined in (|24)) . let us first 
analyze the behavior of /i, i = 0, . . . , n, and /„+i separately. 

Consider Ii, defined in (1^^ . Since, at Di, tk < A — fo'" 
k = i + 1, . . . , n, the functions F{tk \ A — J2j=i+i ^j) ^^'^ continuous, see ([9]). 
The factor F{tn+i \ A — X)j=i+i^i) undergoes a nonzero jump discontinuity 
when point (io, • ■ • , ^n+i) transverses the hyperplane defined as 

n+1 

^ = A, z = 0, . . . , n, (25) 



and is continuous function anywhere else. The result of integration in 
is a continuous function in Z?i, see the proof in Appendix [C] Therefore, at the 
domain Di, each Ii has a discontinuity of a jump type at the hyperplane defined 
in 

Now, consider the continuity of /n+i, expression (|23p. The first term, again, 
is a continuous function in Di , the proof is similar to what is done in Appendix 
[Ul The only discontinuity in the second term at the domain Di is due to the 
factor F{tn+i \ A — X]j=o located at the hyperplane defined as 

n+1 

Y^t.^A, (26) 
j=o 
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while all F{tk \ A — X]^=o ^j)' ^ ~ 0, . . . ,n are continuous functions at this 
domain, see 

According to (HH), the probability density P(t„+i, . . . ,to) can be obtained 
as a sum of all 1^, i — 0, . . . ,n and In+i- Therefore, it inherits all the disconti- 
nuities, contained in li and In+i- So, at the domain Di, the probability density 
P{tn+i, . . . ,to) has nonzero jump discontinuities at the (n + 2) hyperplaneqj 
defined in (j25p and (|26p . and is a continuous function at the rest of the domain. 

4.6 Discontinuities in P{tn+i | tn, • • • ito) 

Conditional probability density P[tn+i \ t„, . . . ,io) can be easily derived from 
([M)) according to the definition It should be outlined, that joint proba- 
bility density P(in, . . . ,<o) is strictly positive for any (n -I- l)-tuple of positive 
values {tm . ■ . ,to) as it can be concluded from ((24| . Moreover, P(i„, . . . , to) is 
continuous at the domain 

n 

J2t^<^■ (27) 

i=0 

Indeed, at the domain ((?7)) . we have also X^r^To^ ^ which means that the 
discontinuities of P{tn, . . . ,to) are located at hyperplanes defined by conditions 
([25)1 and P5|) with (n — 1) substituted instead of n. But those conditions are 
never satisfied due to (|27p . Thus, division of P(i„_|-i, . . . ,to) by strictly positive 
and continuous function P{tn, . . . , to) neither does add new discontinuities, nor 
does it eliminate already found in the P(t„+i, . . . ,to) at the domain Di. 

Therefore, at the domain Di, function P(tn+i \ tn, ■ ■ ■ ,to) contains (n -I- 2) 
jump discontinuities, located at the same positions as in P(t„+i, . . . ,to), equa- 
tions (PS)) and (pS]) , and is a continuous function at the rest of Di . The location 
of discontinuity (1261) depends on to- This dependence cannot be compensated by 
any summands, continuous at hyperplane (1261) . therefore, the whole conditional 
probability density P{tn+i \ t„,...,to) depends on Iq. This means, that the 
condition ([2]) does not hold for any n for the output stream of BN with delayed 
feedback. The Theorem 1 is proven. □ 

5 Particular cases 

In the previous sections, we have proven the impossibility to represent the stream 
of output ISI durations for BN with delayed feedback as a Markov chain of any 
finite order. In particular, output ISI stream is neither a sequence of indepen- 
dent random variables, and therefore is non-renewal, nor it is the first-order 
Markovian process. 

In the course of proving Theorem 1, we have obtained the expression for 
P{tn+i, tn, . . . , to) at the domain X]r=o < A in general case of an arbitrary n, 
see (P^ . This allows to calculate the conditional probability density P(t„+i | 
t„, . . . , to) for X]r=o *i < ^ and n = 0, 1, . . .. 

In this section, we consider two particular cases of P(t„+i | t„, . . . , to) when 
n — and n — 1, namely, the single-ISI conditional probability density P(ti | tp) 
and the double-ISI conditional probability density P(t2 | ti,to) and obtain the 

^Note, that all hyperplanes, defined in I I25I I and I I26I I are difl'erent. 
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expressions for P{ti \ to) and P{t2 \ ti,to) for domain ((TH]), as well as for all 
other possible domains, which were omitted in calculations with arbitrary n. 



5.1 Conditional probability density P(ti \ to) 

In order to derive the exact expression for conditional probability density P{ti \ 
to) for neighbouring ISI durations, we take Steps 1-3, outlined in Section [Ol 
for n = 0. In the case of P{ti \ to), there are only two domains, on which the 
expressions should be obtained separately, namely cases to < A and to > A. 
Performing integration in ((T]), one obtains the following expressions for P{ti, to) 
at these domains: 

Pituto) = Fiti\A)Pito), to>A, 

to 



= I A) j F{to I so)ff(so)dso 





A 



+ / ^^(^1 I So - to)F{to I so)/(so)dso, to < A. (28) 



to 

Expressions can be understood as follows. Since to > A, one can be 
sure that the line has time to get free from impulse during to, therefore at the 
moment of next firing (at the beginning of ti ) the impulse enters the line and 
has time to live equal A. In the case of to < A, see (PSI) . two possibilities arise. 
The first term corresponds to the scenario, when the feedback line discharges 
conveyed impulse within time interval to, and the second one represents the 
case when at the beginning of ti the line still keeps the same impulse as at the 
beginning of to • 

Then, using © and PT|) . one obtains: 

P(ti I to) = F{ti I A), to > A, 

to 

= {F{ti I A) J F(to I so)giso)dso+aF{ti \ A - to)F(to | A) 



A 

+ J F{ti I So - to)F(to I so)g(so)dso), < A. (29) 

to 

It should be outlined, that the output ISI probability density P(to) is strictly 
positive and continuous function at the domain < to < A. Indeed, due to 
(P))- ([T^ . the only discontinuity contained in P(to) is placed at to — A, see 
Figure [21(a). 

It can be shown, that the following normalization conditions take place: 

oo oo 

/ dtiP(ti I to) = 1, and / dtoP(ti,to) = P(ti). 



Using ^ and (f^ . one obtains the positions of discontinuities in P(ti | to): 

ti = A, if to > A, (30) 

ti = A, to + ti = A, if to < A. (31) 
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Figure 3: Conditional probability density P{ti \ to) for t = 10 ms, A = 8 ms, A 
= 400 s^^, A'o — 2, ^0=6 ms (left) and to— 11 ms (right), found numerically by 
means of Monte-Carlo method (the number of firings accounted N = 150 000). 



Obviously, expressions ([5T|) could be obtained directly from (^5]) and by 
substituting n — 0. 

As it can be seen from (j30p and pip , the number of jump discontinuities in 
P{ti 1 ^o) and their positions depend on Iq. Therefore, the conditional prob- 
ability density P{ti \ to) cannot be reduced to output ISI probability density 
P{ti). Therefore, the neighbouring output ISIs of BN with delayed feedback 
are correlated, as expected. 

Examples of P{ti \ to), found for two domains numerically, by means of 
Monte-Carlo method (see Section |6] for details), are placed at Figure [3] 



5.2 Conditional probability density P{t2 \ ti,to) 

In order to derive the exact expression for conditional probability density P{t2 \ 
ti,to) for the successive ISI durations, we take Steps 1-3, outlined in Section H?T1 
for n = 1. In the case of P{t2,ti,to), there are five domains, on which the 
expressions should be obtained separately, namely, the domain 

Di = {{to,ti,t2) I ti+to< A}, 

which was already utilized in Section 21 and the four remaining: 

D2 = {{to,h,t2) \ to>A and ti>A}, 

D3 ^ {ito,ti,t2) \ to<A and ti > A}, 

Di = {{to,ti,t2)\ to>A and ti < A}, 

D5 ^ {ito,ti,t2) \ to<A and A-to<ti<A}, 

Expressions for P{t2 \ ti,to) can be found exactly on each domain: 

P{t2 I h,to) = F{t2 I A), {toMM) e D2, 

^F{t2\A) {to,ti,t2) e D3, 

= F{t2\A~ti), {to,ti,t2) e Di, 
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02468 10 02468 10 

t2, ms t2, ms 

Figure 4: Conditional probability density P{t2 \ ti,io) for t = 10 ms, A = 6 
ms, A — 400 s^^, A'o — 2, ti=8 ms, to=8 ms (left) and ti = 3 ms, = 8 ms 
(right), found numerically by means of Monte-Carlo method {N = 150 000). 



F{t2 I A - h)F{ti\A) f " F{to I sn)5(so)dso 
Jo 



F(t2|A) j F{ti\so - io)F(to|so)5(so)dso + a F{t2\A)F{ti\A - to)F{to\A) 

Jto 

. p^^^ (^(^2 I A - U)F{t, \A) jj F{t„ I so)9{so)dso 



rUs+ti 

+ F{t2 I A) / F{h I So - io)-F(io | so)5(so)dso 



A 

F(t2|so - to - ii)-F'(ii|so - to)-F'(^o|so)5(so)dso 

to+ti 

+ a F{t2\A-to-h)F{h\A-tQ)F{to\A)y {toMM) Di. (32) 

where F(ii,to) = F{h \ A) F{to \ so).9(so)dso + Jt^ F{h \ so - to)F{to \ 
so)/(so)dso, according to ((29|) . 

It is worth to notice, that P(ti , to) is strictly positive and continuous function 
on both Di and D^, see denominators in p2p . Indeed, from ([5(1)) and pil) one 
can see, that P{t i^Iq) may include discontinuities only at the points ti ~ A and 
ti — A — tg. None of these points fall into Di, or D5. 

It can be shown, that the following normalization conditions take place: 

00 00 

/ dt2P{t2 I h,to) = 1, and / dtaP{ta,ti,t2) = P{t2,ti). 


Using © and ([5^ . one derives the positions of jump discontinuities in the 

conditional probability density P(t2 \ ti^to): 

t2^A, {to,ti,t2) e D2UD3, (33) 

h+t2=A, {to,ti,t2) e D^. (34) 

t2 = A, ti+t2=A {to,ti,t2) e D5, (35) 

i2 = A, ti+t2=A, to+ti+t2 = A, {to,ti,t2) e Di. (36) 

Obviously, expression ([55)) could be obtained directly from (ps)) and by 
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t2, ms t2, ms 

Figure 5: Conditional probability density P(t2 \ ti,io) for t = 10 ms, A = 6 
ms, A — 400 s~^, A^o = 2, ti—3 ms, to=3.5 ms (left) and ti — 3 ms, to — 2.5 ms 
(right), found numerically by means of Monte-Carlo method {N = 150 000). 



substituting n — 1. 

As one can see, the number and the position of jump discontinuities in P{t2 \ 
ti,to) depends on to, therefore P{t2 \ ti,to) cannot be reduced to P{t2 \ ti), 
which means that the output stream is not first-order Markovian. 

Examples oi P{t2 \ ti,to), found numerically for different domains, are placed 
at Figures m and [51 

6 Numerical simulation 

In order to check the correctness of obtained analytic expressions, and also 
to investigate whether the output ISIs stream is non-Markovian for inhibitory 
BN with higher thresholds as well as for iVo = 2, numerical simulations were 
performed. A C-I--I- program, containing class, which models the operation 
manner of inhibitory BN with delayed feedback, was developed. Object of this 
class receives the sequence of pseudorandom numbers with Poisson probability 
density to its input. The required sequences were generated by means of utilities 
from the GNU Scientific Librarjl^ with the Mersenne Twister generator as source 
of pseudorandom numbers. 

Program contains function, the time engine, which brings system to the 
moment just before the next input signal, bypassing moments, when neither 
external Poisson impulse, nor impulse from the feedback line comes. So, only 
the essential events are accounted. It allows one to make exact calculations 
faster as compared to the algorithm where time advances gradually by adding 
small time-steps. 

The conditional probability densities, P{ti \ to) and P{t2 \ ti,to), are found 
by counting the number of output ISI of different durations and normalization 
(see Figures [3] - H]) . Obviously, for calculation of conditional distributions only 
those ISIs are selected, which follow one or two ISIs of fixed duration, to for P{ti \ 
to) and {ti, to} for P{t2 \ ti, to). The number and the positions of discontinuities, 
obtained in numerical experiments for inhibitory BN with threshold 2, coincide 
with those predicted analytically in (I5T]) and (155)) - 

^http: / /www. gnu.org/software/gsl/ 
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Figure 6: Conditional probability density P{t2 \ ti,io) for t = 10 ms, A = 6 
ms, A — 1000 s~^, iVo = 4, ti=3 ms, to=3.5 ms (a) and ii = 3 ms, to — 2.5 ms 
(b), found numerically by means of Monte-Carlo method {N — 150 000). 



For A^o > 2, conditional probability densities P{ti \ to) and P{t2 \ ii,io) are 
similar to those, found for No~2. In particular, both the quantity and position 
of discontinuities coincide with those obtained for inhibitory BN with threshold 
2, as expected, compare Figures H] and O 



7 Conclusions and discussion 

Our results reveal the influence of delayed feedback presence on the neuronal 
firing statistics. In the contrast to the cases of BN without feedback [3^ and 
BN with instantaneous feedback [3T , the neighbouring output ISIs of inhibitory 
BN with delayed feedback are mutually correlated. This means that even in the 
simplest possible recurrent network the output ISI stream cannot be treated as 
a renewal one. 

The non-renewalness of experimentally registered spike trains was observed 
for neuronal activity in various CNS areas in mammals [51 |3H1 [33] and fish [51 [7] . 
The simplest stochastic processes which are not renewal are the Markov pro- 
cesses of various order. The order of underlying Markov process was estimated 
in [7] for activity in the weakly electric fish electrosensory system. It was found 
in [7] that for some neural fibers the Markov order should be at list seven, 
which does not exclude that the genuine order is higher, or that the activity is 
non- Markovian . 

Actually, for proving based on experimental data that a stochastic activity 
has Markov order m, one needs increasing amount of data with increasing m. 
If so, it seems impossible to prove experimentally that a stochastic activity is 
non-Markovian. Similarly as it is impossible to prove experimentally that a 
number is irrational. We prove here that the output ISI stream of inhibitory 
BN with delayed feedback is non-Markovian based on complete knowledge of 
the mechanism which generates the output stream. In a sense, to have this 
knowledge is equivalent as to have an unlimited amount of experimental data. 

It is worth to notice, that the activity of excitatory BN with delayed feedback 
is non-Markovian as well [10]. We conclude, that it is namely the delayed 
feedback presence, which results in non-Markovian statistics of neuronal firing. 
One should take this facts into account during analysis of neuronal spike trains 
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obtained from any recurrent network. 



A Proof of Lemma [T] 

In the compound event (t„+i, Sn+i), the time to hve Sn+i always gets its value 
before than the tn+i does. The value of Sn+i can be determined unambiguously 
from the (t„,s„) value (See Sections [2.21 and : 

The only two factors, which determine the next ISI duration, are (i) 

the value of Sn+i, and (ii) the behavior of the input Poisson stream under 
the condition (t„, s^; . . . ;to, sq) after the moment 0, when the starts. The 
Sn+i value does not depend on (i„_i, Sn-i] • ■ • ; sq), see above. As regards the 
input Poisson stream, condition s„; . . . ; to, so) imposes certain constraints 
on its behavior before the 9. Namely, if ti ^ Si for some < i < n, than one 
can conclude that an input impulse was obtained just at the end of ti . In the 
opposite situation, when ti — Si, one can conclude that in the course of ti exactly 
one impulse was obtained from the Poisson stream. But what do we need in the 
definition of the P(t„+i, s„+i | s„; . . .;to, Sq), it is the conditional probability 
to obtain input impulses at definite moments after the 9. For a Poisson stream 
this conditional probability does not depend on conditions before the 9. For 
example, conditional probability to obtain the first after 9 impulse at 6* + t 
equals e~^^Xdt, whatever conditions are imposed on the stream before the 9. 
This proves □ 



B Finding integrals Ii for P{tn+i, ■ ■ ■ ,to) 

Domain of sq values covered by Ii, i = 0, . . . , n, corresponds to the scenario, 
when impulse, which was in the feedback line at the beginning of interval to 
(with time to live sq), will reach BN during interval ti, see Figure [T] In this 
process, after each firing, which starts ISI tk, k < i, the time to live of the 
impulse in the feedback line is decreased exactly by tk-i- This means, that 
variables of integration {sq, . . . , Sn+i}, above, are not actually independent, but 
must satisfy the following relations: 

fe-i 

Sk = so-^tj, = (37) 

j=o 

which are also ensured by 5-function in the bottom line of (fT7)) . Next to Si time 
to live must be equal A: 

S^+l = A, (38) 

and this is ensured by (5- function in the top line of (jl7|) . 

The next to s^+i times to live again are decreased by corresponding ISI with 
each triggering. Due to p^ . this brings about another set of relations: 

k-l 

Sfe = A- ^ tj, k^i + 2,...,n + l, (39) 
j=i+i 
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to tl ti—i ti ti+i t„ time axis 

Figure 7: Illustration of relations between (fo, . . . , t„) and (sq, • • • , Sn+i) con- 
tributing to the li'. sq G ] X]j=o ^J' Sj=o ^j] ' 5I]J=o ^^'^ time to live 
decreases steadily with every output firing for A; = 0, i — 1 until it becomes 
that Si <ti. Then, during the time interval ti the line discharges its impulse to 
BN input, and at the beginning of U+i starts to convey the new one with time to 
live Sj+i = A. After that, times to live Sk are again decreased by corresponding 
tk with each firing, k = i + l, n. 
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which are again ensured by J-function in the bottom hne of PT|) . Relations ([57]) , 
([55)1 and ([55)1 together with limits of integration over so in ([T51) ensure that at 
Di the following inequalities hold: 

Sk >tk, /c = 0, . . . ,i - 1, 

Sk>tk, k = i + l,...,n. (40) 

Inequalities PO)) allow one to decide correctly which part of rhs of (IT7)) should 
replace each transition probability P{tk, Sk \ tk~i, Sfe-i) in (jl9p . and perform all 
but one integration. This gives: 

li^ / dso / dsi • . . . • / ds„+i 



i fe-1 

i^(io I so)/(so) n I - ^0 + 5Z 

fc=l j=0 

n+1 fe-1 

X F{U+i I s,+i) (5(s,+i - A) F{tk I Sfc)(5(sfc - A + ^ t^) 

J^' t ■ 

n+1 fc-1 ^J-" ' i 

= Y[ F{tk\A- h) / X{F{tk\sa~Y,tj)g{sQ)Aso, 

k=i+l j=2+l ^1 fe=0 j=0 

i = 0,1,2,..., n. (41) 

The last expression might be obtained as well by means of consecutive sub- 
stitution of either top, or bottom line of (jl7D into (TT^ . without previously 
discovering ^7} - (|40)) . 

Finally, integral /„+i corresponds to the case, when at the beginning of 
interval tn+i, the line still keeps the same impulse as at the beginning of io- 
Therefore, comprises the rest of scenarios contributing to the value of 

P(tn+i, ■ . ■ jtg) in ([5]). Proceeding as in the preceding terms, the contribution 
Ji+i reads: 

^A ^A 
In+i = dso / dsi . . . / ds„+i 

n+1 fe-1 

-F(io I so)/(so) n ^(^fe I ^k)5{sk - So + ^ij) 

fc=l j=0 

ri+1 fe-1 



Y[F{tk\so-J2tj)f{s„)dso 

^n. f k=0 j=0 

n+1 fe-1 n+1 fe-1 

n ^i^k I so - ^ tj)g{sQ)dso + a n ^(^^ ' ^ " H (^2) 

J n n 7, r\ a n 



Y-n fc=o i=o fe=o i=o 
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C Continuity of integral factor in (1221) 

Continuity in Di of the integral factor 

^y*^' ^ k-1 

^ X!^-J')f(so)dso, i = 0,l,...,n, (43) 

1, k=0 3=0 

in the expression (j22D can be proven after mathematical simplification. First, 
notice that due to integration domain the following inequalities take place 

k-l i-1 

so-^tj>tk, fc = 0, 1, . . . ,i - 1, So ^y^Jj < U, 

which together with © allows to replace with the following 

i-1 ^^'f^' 

n (A'ifee-^*'=) / (l + Asi)e-^^ip"(i,-si)5(so)dso, 

where Si = Sq — X]}=o ■ The continuity of the last expression is determined 
by the continuity of its second factor, since the first one is continuous in M"+^. 
The second factor can be replaced with 

^(1 + Xso)e-~^'^^P"{U - so)g (^so + fc^J j ^sq. (44) 

after changing the variable of integration. For further simplification of the last 
expression use ([1]), (IT^ and pH)) . which gives instead of pi| 



= e-^*H, j{l + Xso)X^g ( so + ^ j dsQ- (45) 





ti 



i-1 



-e-^*' J{1 + Aso)A2so g + ^ tjj dso. (46) 

The required continuity of (j43p is determined by the continuity of integral factors 
in (|^5|) and Now, take into account the explicit expression for g{s), which 
is found in pMi Eq. (15)]. For our purposes it is enough to know that g{s) — 
A + Be'^'^'^ , where A and B are constants. Taking this into account, the integral 
factor in P5)) can be replaced with 

ti ti 

A [{l + \so)\^dso + Be^^^r=ot3 [ {1 + \so)\^e^^'»dso, 



which makes its continuity self-evident. The same is for integral factor in (1461) . 
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